In [1]:
%pylab inline
import matplotlib.pyplot as plt

import collections

Populating the interactive namespace from numpy and matplotlib

Band Structure Calculations by Pseudopotential Methods

for Si, diamond structure

Peter Y. Yu, Manuel Cardona: Fundamentals of Semiconductors Physics and Materials Properties see here

Further reading:

Rolf Enderlein & Norman J Horing: Fundamentals of Semiconductor Physics and Devices see here

In [2]:
# Az abra kimentesehez az alabbiakat a plt.show()  ele kell tenni!!! 

#savefig('fig_rainbow_p2_1ray.pdf');  # Abra kimentese
#savefig('fig_rainbow_p2_1ray.eps');  # Abra kimentese

# Abra es fontmeretek
xfig_meret= 7  #    12 nagy abrahoz
yfig_meret= 8    #   12 nagy abrahoz
xyticks_meret= 15  #  20 nagy abrahoz
xylabel_meret= 20  #  30 nagy abrahoz
legend_meret= 21   #  30 nagy abrahoz
In [3]:
def Ham_op(kv, Gvec, a, Vpseudo_form_s, Vpseudo_form_a):
    
    # pesudo pot. 
    # energia in units of Rydberg=13.6 eV = hbar^2/(2m a_0^2), 
    #where a_0 is the Bohr radius
    ##  Gvec ---> reciprocal vectors
    ## a ---> lattice constant
    # Vpseudo_form_s, Vpseudo_form_a ---> Local Pseudopotential Form Factors: 
    # see: J. R. Chelikowsky and M. L. Cohen, Phys Rev. B 10, 12 (1974). 
    # and Peter Y. Yu, Manuel Cardona: Fundamentals of Semiconductors 
    # Physics and Materials Properties
    
    # Hsize = the size of the Ham matrix
    
    aB= 0.053  #  Bohr radius (in units of nm),  a_0=5.29×10−11 m
    #a = 0.5431 #  lattice constant for Si (in unit of nm),   5.431 A
    
    fact = (2*pi/(a/aB))**2  #  to get the kinetic energy in units Rydberg
    
    s= 1/8*array([1,1,1])
    Hsize = len(Gvec)  # initialization of Ham
    Ham = matrix(zeros((Hsize,Hsize),complex))
    for i in range(Hsize):
        for j in range(Hsize):
            Gij=  Gvec[i] - Gvec[j]
            tmp = dot(Gij,Gij)
            
            if tmp in [0,3,4,8,11]:
                Vs = Vpseudo_form_s[tmp]
                Va = Vpseudo_form_a[tmp]
            else:
                Vs=0
                Va=0
            
            g_dot_s = dot(2*pi*Gij,s)
            Ham[i,j] = Vs*cos(g_dot_s)-1j*Va*sin(g_dot_s)
            #Ham[i,j] = 0   ##  free electron
    
        # a Ham diagonalis elemei = a kinetic energy
        g = Gvec[i]
        Ham[i,i] = fact* dot(kv-g,kv-g) 
      
        # Ham is in units of Rydberg
        
    return Ham

Local Pseudopotential Form Factors:

see: J. R. Chelikowsky and M. L. Cohen, Phys. Rev. B 10, 5095 (1974). https://journals.aps.org/prb/abstract/10.1103/PhysRevB.10.5095

and Cohen & Chelikowsky: Electronic Structure and Optical Properties of Semiconductors

http://www.springer.com/la/book/9783642970801

Form Factor (Ry) Si Ge
$V_3$ -0.2241 -0.2768
$V_8$ 0.0551 0.0582
$V_{11}$ 0.0724 0.0152
In [4]:
##
a_Si = 0.5431 #  lattice constant for Si (in unit of nm),   5.431 A
a_Ge = 0.5657 #  lattice constant for Si (in unit of nm),   5.431 A

### pseudo potential form factors for Si

#Vs_Si = {0:0, 3:-0.211, 4:0, 8:0.04, 11:0.08}
#Va_Si = {0:0, 3:0, 4:0, 8:0, 11:0}  

Vs_Si = {0:0, 3:-0.2241, 4:0, 8:0.0551, 11:0.0724}
Va_Si = {0:0, 3:0, 4:0, 8:0, 11:0}  

Vs_Ge = {0:0, 3:-0.2768, 4:0, 8:0.0582, 11:0.0152}
Va_Ge = {0:0, 3:0, 4:0, 8:0, 11:0}  

#Vs_Si[3],Va_Si[8],Vs_Si[0],Va_Si[0]
In [5]:
li = ['a', 'b', 'new', 'mpilgrim', 'z', 'example', 'new', 'two', 'elements']
if('example' in li):
    print("fifi")

fifi

Some results here:

path: specpoints = array([Lp,Gp,Xp,Up,Gp])

Compare with figure in Peter Y. Yu, Manuel Cardona: Fundamentals of Semiconductors Physics and Materials Properties, Fig. 2.10, page 73 (in pdf version)

-----------------------------

path: array([Gp,Xp,Up,Gp,Lp,Kp,Wp])

file: Lecture2_band_struct.pdf

-----------------------------

Aaron Danner (2011),

path: array([Lp,Gp,Xp,Wp, Kp,Gp])

https://www.ece.nus.edu.sg/stfpage/eleadj/pseudopotential.htm

In [6]:
# Diamond lattice, 'a' is the lattice constant of the fcc cube

## unit cell vectors of the reciprocal vectors (in units of 2 pi/a):

b1 = array([-1,1,1])
b2 = array([1,-1,1])
b3 = array([1,1,-1])

# creating the reciprocal vectors (in 'array' form): 
Gnum = 3  #### (2*Gnum+1)**3, number of Gvec
indx=array(range(-Gnum,Gnum+1))
Gvec = []
for n1 in indx:
    for n2 in indx:
        for n3 in indx:
            Gvec.append(n1*b1+n2*b2+n3*b3)

# special points in the Brillouin zone:
Gp = array([0,0,0])
Xp = array([0,1,0])
Wp = array([1/2,1,0])
Kp = array([3/4,3/4,0])
Lp = array([1/2,1/2,1/2])
Up = array([1/4,1,1/4])

## FIGYELEM: fcc racs eseten az irodalomban az U,K jeloles azt jelenti, hogy 
## U-K vonal menten NEM rajzoljuk ki a savot!!!!!

###  kiiratashoz, az x-tengelyen:
specNev_map = {tuple(Gp):"$\Gamma$",tuple(Xp):"$X$",tuple(Wp):"$W$",
               tuple(Lp):"$L$",tuple(Kp):"$K$",tuple(Up):"$U$"}
In [7]:
## Si with diamond lattice structure

Npoints = 35;  ## hany pontbol keszul a gorbe
ymin=-14
ymax = 6.2;

figsize(xfig_meret,yfig_meret)

#### path through special points in the Brilloin zone
#specpoints = array([Gp,Xp,Wp,Lp,Gp,Kp,Xp]) 
#specpoints = array([Lp,Gp,Xp,Up,Kp,Gp])    
     
#specpoints = array([Gp,Xp,Up,Gp,Lp,Kp,Wp]) 
#specpoints = array([Lp,Gp,Xp,Up,Gp])  
specpoints = array([Lp,Gp,Xp,Wp, Kp,Gp])   # Aaron Danner, 2011  

klist = []
specpos = [0]
kk = [0]
tmp = 0
for i in (range(len(specpoints)-1)):
    sl = specpoints[i+1]-specpoints[i]  
    sl_hossz = sqrt(dot(sl,sl))
    dk = sl/Npoints
    for num in arange(0,Npoints):
        # k vector along the line given by speclines:
        klist.append(specpoints[i] + num*dk)
        # length of the k vector and shifted: 
        tmp +=  sqrt(dot(dk,dk)) 
        kk.append(tmp) 
    specpos.append(tmp)

klist.append(klist[-1]+dk)   

enlist = []
for kv in klist:
    enlist.append(eigvalsh(Ham_op(kv, Gvec, a_Si, Vs_Si, Va_Si)))
    
E_shift = 10.52211319  ## in units of eV
# az energiskala ugy van eltolva, hogy a G pontban a legfelso valencia sav erteke = 0

enlist= 13.6*array(enlist) - E_shift  ###  in unit of eV 

eigszam= len(enlist[0,:])
print("The # of eigenvalues = ",eigszam,"= # of G vectors =",len(Gvec))
      
# drawing the figure

n_eigval=len(enlist[0,:])
n_eigval=10
for ii in range(0,n_eigval):
    plot(kk,enlist[:,ii],color='r')

i = 0
for xc in specpos:
    axvline(x=xc,color='b')
    text(xc,-15,specNev_map[tuple(specpoints[i])],fontsize=xylabel_meret)
    i=i+1
    
ylabel(r'$E(k)$ (eV)',fontsize=xylabel_meret)

#  kikapcsoljuk a xticks: 
xticks([], [])

xlim(0,specpos[-1])
ylim(ymin,ymax)

yticks(arange(-14,5,2))

cim = "Si with Local Pseudopotential"
title(cim,fontsize=xylabel_meret)

grid();

#savefig('Fig_pseudo_LGXWKG.png');  # Abra kimentese

The # of eigenvalues =  343 = # of G vectors = 343

Gap = $E_3(X)$-$E_2(\Gamma)$

In [8]:
E1= 13.6*eigvalsh(Ham_op(Gp, Gvec, a_Si, Vs_Si, Va_Si))[1]
E2= 13.6*eigvalsh(Ham_op(Xp, Gvec, a_Si, Vs_Si, Va_Si))[4]
Gap = E2-E1
print("E(G) = %5.2f, E(X) = %5.2f, GAP (eV) = %5.2f " % (E1, E2, Gap))

E(G) = 10.27, E(X) = 11.43, GAP (eV) =  1.17
In [9]:
13.6*eigvalsh(Ham_op(Gp, Gvec, a_Si, Vs_Si, Va_Si))[:11]
Out[9]:
array([ -2.32872833,  10.26657191,  10.26657191,  10.26661498,
        13.62936999,  13.62936999,  13.62941114,  14.40551073,
        18.03363896,  18.03363896,  18.63782958])

U,K points together !!!!!!!!

In [10]:
## Si with diamond lattice structure

Npoints = 40;  ## hany pontbol keszul a gorbe
ymin=-14
ymax = 6.2;

# creating the reciprocal vectors (in 'array' form): 
Gnum = 2  #### (2*Gnum+1)**3, number of Gvec
indx=array(range(-Gnum,Gnum+1))
Gvec = []
for n1 in indx:
    for n2 in indx:
        for n3 in indx:
            Gvec.append(n1*b1+n2*b2+n3*b3)
            
figsize(xfig_meret,yfig_meret)

#### path through special points in the Brilloin zone
#specpoints = array([Gp,Xp,Wp,Lp,Gp,Kp,Xp]) 
#specpoints = array([Lp,Gp,Xp,Up,Kp,Gp])    
#specpoints = array([Gp,Xp,Up,Gp,Lp,Kp,Wp]) 
#specpoints = array([Lp,Gp,Xp,Wp, Kp,Gp])   # Aaron Danner, 2011  
specpoints = array([Lp,Gp,Xp,Up,Kp,Gp])

## FIGYELEM: fcc racs eseten az irodalomban az U,K jeloles azt jelenti, hogy 
## U-K vonal menten NEM rajzoljuk ki a savot!!!!!

###  kiiratashoz, az x-tengelyen:
specNev_map = {tuple(Gp):"$\Gamma$",tuple(Xp):"$X$",tuple(Wp):"$W$",
               tuple(Lp):"$L$",tuple(Up):"$U,K$",tuple(Kp):"$\Gamma$"}



klist = []
specpos = [0]
kk = [0]
tmp = 0

for i in (range(len(specpoints)-1)):
    if not i==3:
        sl = specpoints[i+1]-specpoints[i]  
        sl_hossz = sqrt(dot(sl,sl))
        dk = sl/Npoints
        for num in arange(0,Npoints):
           # k vector along the line given by speclines:
             klist.append(specpoints[i] + num*dk)
           # length of the k vector and shifted: 
             tmp +=  sqrt(dot(dk,dk)) 
             kk.append(tmp) 
        specpos.append(tmp)

klist.append(klist[-1]+dk)   

enlist = []
for kv in klist:
    enlist.append(eigvalsh(Ham_op(kv, Gvec, a_Si, Vs_Si, Va_Si)))
    
E_shift = 10.52211319  ## in units of eV
# az energiskala ugy van eltolva, hogy a G pontban a legfelso valencia sav erteke = 0

enlist= 13.6*array(enlist) - E_shift  ###  in unit of eV 

eigszam= len(enlist[0,:])
print("The # of eigenvalues = ",eigszam,"= # of G vectors =",len(Gvec))
      
# drawing the figure

n_eigval=len(enlist[0,:])
n_eigval=10
for ii in range(0,n_eigval):
    plot(kk,enlist[:,ii],color='r')

i = 0
for xc in specpos:
    axvline(x=xc,color='b')
    text(xc,-15,specNev_map[tuple(specpoints[i])],fontsize=xylabel_meret)
    i=i+1
    
ylabel(r'$E(k)$ (eV)',fontsize=xylabel_meret)

#  kikapcsoljuk a xticks: 
xticks([], [])

xlim(0,specpos[-1])
ylim(ymin,ymax)

yticks(arange(-14,5,2))

cim = "Si with Local Pseudopotential"
title(cim,fontsize=xylabel_meret)

grid();

#savefig('Fig_pseudo_LGXU-KG.png');  # Abra kimentese

The # of eigenvalues =  125 = # of G vectors = 125
In [ ]: